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ABSTRACT 

The observed wide eccentricity distribution of extrasolar giant 
planets is thought to be the result of dynamical instabilities and 
gravitational scattering among planets. Previously, it has been 
assumed that the orbits in giant planet systems become gravita- 
tionally unstable after the gas nebula dispersal. It was not well 
understood, however, how these unstable conditions were estab- 
lished in the first place. 

In this work we numerically simulate the evolution of sys- 
tems of three planets as the planets sequentially grow to Jupiter's 
mass, and dynamically interact among themselves and with the 
gas disk. We use the hydro-dynamical code FARGO that we mod- 
ified by implementing the iV-body integrator SyMBA. The new 
code can handle close encounters and collisions between planets. 
To test their stability, the planetary systems were followed with 
SyMBA for up to 10 8 yr after the gas disk dispersal. 

We find that dynamics of the growing planets is complex, be- 
cause migration and resonances raise their orbital eccentricities, 
and cause dynamical instabilities when gas is still around. If the 
dynamical instabilities occur early, planets can be removed by 
collisions and ejections, and the system rearranges into a new, 
more stable configuration. In this case, the planetary systems 
emerging from the gas disks are expected to be stable, and would 
need to be destabilized by other means (low-mass planets, plan- 
etesimal disks, etc.). Alternatively, for the giant planet system 
to be intrinsically unstable upon the gas disk dispersal, a special 
timing would be required with the growth of (at least some of) 
the giant planets having to occur near the end of the gas disk 
lifetime. 

Key words: Giant planets dynamics, hydro-codes, N-body sim- 
ulations. 



1 INTRODUCTION 

The eccentricity distribution of the extrasolar gi- 
ant planets is wide with orbits commonly having 
e > 0.3. Such a wide distribution was unexpected 
based on our anticipation from the Solar System 
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planets. Different mechanisms have been pro- 
posed to explain the high eccentricities values: 
trapping of planetary pairs in mean motion res- 



teractions with the gas. According to the Pollack 



onances (Lee and Peale (2002 i), the Kozai cycles 



in binary systems ( Mazeh et al. ( 1997 I ; |Holman 
et Etl.j ( |1997"1 )), stellar jets ( |Namouni| ( |2005| >), and 
gravitational scattering during global dynamical 



instabilities ( Weidenshilling and Marzari ( 1996 I; 



Rasio and Ford (1996); Lin and Ida (19971) 



et al. (20081; Raymond et al. (20091; Juric and 



This last mechanism has been investigated 
with extensive iV-body simulations ( |Chatterjee 



Tremaine ( 2008 1). All these studies assumed that 
the planetary systems emerging from the gas 
disks are intrinsically unstable, and the gravi- 
tational interactions among planets cause insta- 
bilities after the gas disk dispersal. The subse- 
quent scattering encounters between planets lead 
to large orbital eccentricities, just as needed to 
explain the observations. 

A hydro-dynamical code has been recently 
used to study the planetary system instabilities 



in a low-density rapidly-dispersing disk ( Moeckel 



|and A rmitagc (2012)). The initial orbits of the 
planets were assumed to be circular and close 
enough to each other to be unstable, without 
any mutual resonance relationship. It was found 
that the disk can stabilize some of the plan- 
etary systems by driving them into resonance 
rapidly. However, the systems that became un- 
stable ended up behaving as in the gas-free sim- 
ulations. 

The investigations discussed above, includ- 



ing Moeckel and Armitage (20121, adopt similar 



choices of initial conditions with unstable and 
sometimes overlapping planetary orbits. In real- 
ity, the initial conditions of these studies should 
be informed from the previous stages of planet- 
disk interactions when the damping effects of gas 
were important. 

The orbital dynamics of giant planets in a 
massive gas disk has been studied with a hydro- 
dynamical code by Marzari et al. (20101. They 



started their simulations with the fully-formed 
giant planets and ignored the previous stage dur- 
ing which the giant planets had grown by gas 
accretion onto their cores. 

Here we report the result of the first effort 
to investigate the dynamical evolution of plan- 
ets from their growth phase from cores to after- 
math of the gas disk dispersal. Our simulation 
set up is different from those of previous works 
and is defined according to the following ratio- 
nale. Planetary cores are expected to form by 
oligarchic growth ( jKokubo fc Ida| ( |1998[ )) with 
orbital separations of about 10 mutual Hill radii. 
However, during and after their formation they 
migrate in the disc due to their gravitational in- 



et al. model (Pollack et al. (1996)) the cores can 
spend a few millions years in the disc before ac- 
creting gas in a runaway fashion and become gi- 
ant planets. In this time they can substantially 
modify their orbits and reach a new equilibrium 
configuration. While it was thought in the past 
that cores continuously migrate toward the cen- 



tral star (Ward (19971), it is now known that in 



a disc with realistic heat diffusion they migrate 
towards an orbital radius where migration is can- 



et al. 



celled ( Paardekooper & Mellema (20061; Kley 



(20091; Lyra, Paardekooper, & Mac Low 



( 2010 1) . This no-migration radius acts as a planet 
trap. If multiple cores are present they are ex- 
pected to reach a resonant non migrating config- 



uration near the trap (Morbidelli et al. ( 2008 1) . 



In this configuration the cores can be much closer 
to each other than their initial 10 Hill radii sepa- 
ration which may lead to very strong instabilities 
when the planets grow to Jupiter mass. 

With this kind of dynamical evolution in 
mind, in our hydrodynamical simulations we set 
up a planet trap and let a system of three em- 
bryos of 10 Earth masses to evolve until they 
reach a stable resonant configuration. Then, we 
track the evolution of the systems as each of the 
three planets grows in sequence to one Jupiter 
mass. Finally, we slowly remove the gas from 
the disc and follow the evolution of the sys- 
tems up to 10 s years after the gas dispersal. 



We use the hydro-dynamical code FARGO ( Mas- 
set (2000 1 ) modified in Morbidelli and Nesvorny 



(20121 b y implementing the iV-body integrator 
SyMBA jDuncan et"aL](|l998[))P1to handle close 



encounters and mutual collisions between plan- 
ets. 

The paper is organized as follows. Section 
2 explains the set-up of our numerical simula- 
tions. The dynamics of growing planets is de- 
scribed in Section 3. We then discuss the effects 
of the gas disk dispersal and the subsequent stage 
of purely TV-body interaction of the remaining 
planets (Section 4). Conclusions are given in Sec- 
tion 5. 



2 SETUP OF NUMERICAL 
SIMULATIONS 

2.1 Disk parameters 

We used the hydro-dynamical two-dimensional 



code FARGO (Masset (2000)), in which the orig- 



1 We use swif t_symba7 that is capable of correctly 
handling the closely-packed planetary systems ( |Lev-| 
|ison et al.| ( |2011[ |). 
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inal JV-body Runge-Kutta integrator was re- 



placed ( Morbidelli and Nesvorny (2012l) with 
the symplectic integrator SyMBA (Duncan et 



al. (19981). The SyMBA code was specifically 



designed to handle close encounters and mu- 
tual collisions between planets. As the hydro- 
dynamical simulations are CPU expensive, we 
were not able to run many simulations to fully 
explore parameter space. Instead, we considered 
a few cases that illustrate different aspects of the 
problem. 

Two different disks were considered (de- 
noted by A and B in the following) with the 
initial mass Ma = 0.Q09.M© in case A and 
Mb = O.OI8M0 in case B, where Mq is the mass 
of the Sun. In each case we performed several 
simulations that differed in the prescription for 
the growth of the planets (see below). 

We use units such that G = 1 and Mq = 1. 
The orbital period of a planet with semi-major 
axis a = 1 is therefore T = 2n. We normalize 
the time t by T in the following, so that t cor- 
responds to the number of orbits at a — 1, or 
years. The disk's kinematic viscosity coefficient 
is set to be v — 10~ 5 in these units. The initial 
surface density profile scales with the distance r 
from the star as r'- -1 '' 2 - 1 . 

Our computational domain consists of an 
annulus of the protoplanetary disk extending 
from r min to r max . Different disk extensions have 
been used in different cases. In some cases (spec- 
ified below), the disk had to be extended during 
the simulations when the migration caused the 
innermost planet to approach the inner bound- 
ary. To start with we used 

4.5, and a grid of N r = 660 linearly spaced ra- 
dial cells and N s = 700 azimuthal cells. 

The width of planet's horseshoe region is 
given, in the isothermal disk approximation 



( Masset et al. [2006b l), by 



(m/Me) 

(H/r) a - W 

For example, for a planet of mass tti/Mq = 
3 x 10 -5 (i.e. 10 Earth masses) and disk aspect- 
ratio H/r — 0.05 we get x — 0.0245a, where a 
is the semi-major axis (see Section 2.2). The ra- 
dial resolution of 0.006 allows us to resolve the 
horseshoe region by at least 4 cells for any a ^ 1 . 

The planetary contribution to the potential 
«J> acting on the disc is smoothed according to: 



$ = -- 



Gm 



(2) 



where d denotes the distance of a disc element 
to the planet and e is the smoothing-length. In 
our simulations we used e = 0.5Rh, where Rh 
denotes the Hill radius. 



2.2 Initial orbits 

Three 10-Earth-mass planetary cores were 
placed into the disk and were initially evolved 
till they reached a stable configuration in a res- 



onance. We used a planet trap (Masset (20021; 



Masset et al.| ( |2006[ ); |Morbidelli et al.| ( |2008[ )) to 
halt the orbital migration of the innermost core. 

The planet trap was set as a steep and lo- 
cally positive surface-density gradient in the disk 
inside the initial orbital location of the innermost 
core. It allows the system of three cores to ac- 
quire stable, separated and non-migrating orbits. 
The planet trap is convenient way to mimic the 
situation in real radiative disks where the non- 
isothermal effects can change the direction of 



the type-I migration ( Paardekooper & Mellema 



(2006); Kley et al. (2009)). The migration in the 
inner part of a real radiative disk can be di- 
rected outward, while it remains directed inward 
in the outer disk. This establishes the existence 
of a critical radius where migration vanishes. The 
planetary cores migrating inwards will be col- 
lected near this radius as in the case of a planet 



trap (Lyra, Paardekooper, & Mac Low (2010 1). 

The planet trap location was set at a = 3 in 
case A and a — 2 in case B. The local and posi- 
tive surface-density gradient required to form the 
planet trap was created by imposing a transition 
in the viscosity from 4v to v over Ar = 1 around 
the trap location (Masset et al. (2006 1 ) The ini- 



tial orbits of the three cores were chosen near 
the 5:4 resonant chain in case A (semimajor axes 
ai = 3.07, a,2 = 3.62 and 0,3 = 4.20), and near 
the 3:2 resonant chain in case B (semimajor axes 
Oi = 2.1, a 2 = 2.77 and a 3 = 3.66). The initial 
eccentricities were set to zero. 

In a first step, we followed the evolution of 
disk and cores, and waited till the cores arranged 
themselves in a stable resonant configuration. In 
case A, the cores 3 and 2 ended up in the 6:5 res- 
onance, and cores 2 and 1 in the 7:6 resonance. 
In case B, two cores reached a coorbital configu- 
ration (1:1 resonance) near the planet trap, and 
the third one ended up in the 6:5 resonance with 
the other two. The eccentricities remained small 
at this stage, e ~ 0.01-0.02, due to the strong 
damping of gas, and the orbits remained nearly 
coplanar. 



2.3 Mass growth 

Once the resonant configuration was achieved, 
the mass of each core was increased from the 
initial value (m(0) = 3 x 10 _5 Mq) to one Jupiter 
mass (m.j — 10~ 3 Af©) as follows: 



m(t) = m(0) + (m.j— m(0)) sin 2 (- 



*(0)), 



At 



,(3) 
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where t(0) was the time when the growth started, 
and At was the growth time interval. Different 
values of t(0) were chosen for different planets, 
so that they grew in sequence. Sometimes, we let 
the innermost core grow first with the other two 
growing later. In other cases we opted for grow- 
ing the middle or outer core first (see Section 3). 

We didn't consider gas accretion within the 
Roche lobe. This is a delicate point which is not 
well understood yet and which goes beyond the 
purpose of the present work. 

The criterion for collision is that the dis- 
tance between planets becomes equal or lower 
to one Jupiter radius. 

The timescales for the planetary growth that 
we adopt ( At = 10 3 , 3 x 10 3 , and 4 x 10 3 , see table 
3.1 1 are too short to be realistic. These timescales 
are dictated by our current CPU power (using 
the parallel FARGO code and 30 CPUs we com- 
pute 10 to 100 orbits in 1 hour, depending on 
the disk parameters). Slower growth rates will 
be investigated in the future. 

When a planet grows to Jupiter mass, it is 
expected that the planet trap should become in- 
effective and the planet should start migrating 
inward. We find this behavior in our simulations. 



3 DYNAMICS IN THE GAS DISK 

Here we discuss the orbital dynamics of grow- 
ing planets in the full gas disk. Each system 
is evolved over 20000 to 65000 years after the 
growth of the last planet. To identify the dif- 
ferent settings of our simulations we label them 
Ai t j t k for case A and -Bi,j,fc for case B, where 
indices i,j,k £ [1,3] indicate the growth order 
of the three cores (initially ai ^ a 2 < 0,3). The 
characteristics of each simulation are reported in 
Table O 



3.1 Case A 2 ,i, 3 

Figure [T] shows the results of simulation ^2,1,3 • 
We find that when core #2 grows to Jupiter mass 
it starts migrating inward and scatters the other 
two cores at a > 3. We remark that the initial 
location of the trap is at a — 3, but the gap 
opened by core #2 has shifted the trap at a — 4. 
Therefore, cores #1 and #3 remain trapped at 
quasi constant semi-major axis, respectively at 
a ~ 4 and a ~ 5.6 till core #1 grows. As the 
fully grown planet #1 starts migrating inwards 
at t = 8000, it opens a gap and shifts the trap 
location at a ~ 5.5. Core #3 is scattered out and 
remains near the new trap position at a ~ 5.5 till 
£3(0) = 12000. Core #3 then starts growing and 
inward migration begins. 



Once all three planets reach Jupiter mass 
they migrate inward and evolve into mutual res- 
onances. Their orbital eccentricities rapidly grow 
to large values by resonant interactions. The sys- 
tem becomes highly unstable. The gas density 
distribution is strongly perturbed at this point 
(Fig. [2J leading to complex gas-planet interac- 
tions. Finally, at t = 20000 years, the inner 
planet is lost by plunging into the star. 

Our criterion for collision with the star is 
that the pericenter of the planet's orbit becomes 
smaller than 0.01. The tidal effects are ignored 
in our simulation. In reality, however, the tidal 
effects should start to be dominant for small peri- 
centers, potentially leading to the circularization 
of the orbit in the hot Jupiter region (IBeauge 



& Nesvorny (20121). This effect would result in 
decoupling the planet from the other two. So, in 
any case, the inner planet's influence on the other 
two is suppressed. The subsequent evolution of 
the two-planet system leads to stable orbits with 
moderate eccentricities. 



3.2 Case A 3 , li2 

In our first simulation core #3 was grown on 
At = 10 3 . One of the two remaining cores was 
ejected from the system during the growth of 
core #3. We have therefore reconsidered this 
case with a longer p hase of growth of core #3, 
At 3 = 4 x 10 3 (table 3.1 1. In this case, the tran- 
sition is less violent, core #2 is scattered out 
and its semi-major axis remains stable at the 
new trap location (after gap opening by planet 
#3), i.e. at a ~ 4.3. Core #1 migrates inward 
at the same rate as the fully grown planet #3 
(Fig. [3). The growth of core #1 then leads to a 
phase when the other two planets are scattered 
outward. The subsequent dynamics is complex 
with episodes of outward migration, and a rapid 
increase of the eccentricities after the growth of 
planet #2. As in case ^2,1,3, the system becomes 
highly unstable. Planet #1 is then ejected from 
the system. These early ejections could be re- 



lated to free-floating planets (Sumi et al. (2011 1). 

The orbits of the remaining two planets are 
chaotic till t = 52594, when the two planets 
merge QThe remaining giant planets migrates in- 
ward and converges to a circular orbit. 



2 Note that merging events may happen too often in 
our simulations due to the coplanar approximation 



of the system. 
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Simulation [r mirl ,r max \ 


At 


(N r ,N s ) 





A 2 ,l,3 


[0.5,6.5] 


10 3 


(660,700) 


■^3,1,2 


[0.5,9.5] 


Ai 3 = 4 x 10 3 


(660,700) 






Aii, 2 = 10 3 




-^3,2,1 


[1,9.5] 


3 x 10 3 


(660,700) 




[0.5,9.5] t > 20000 




(700,700) 


Si, 2, 3 


[0.5,4.5] 


10 3 


(660,700) 




[0.1,4.5] t > 8500 




(720,700) 


^3,1,2 


[0.5,4.5] 


10 3 


(660,700) 




[0.1,4.5] t > 8500 




(720,700) 




[0.3,4.5] t > 30000 




(690,700) 



Table 1. Simulations parameters. 



3.3 Case A 3t2 ,i 

In case ^3,2,1, all the three cores grow on At = 
3 x 10 3 (Fig. [4J). The growth phase leads to the 
inward migration of the growing planet #3. The 
other two cores are scattered outward. The sub- 
sequent phase is quite different with respect to 
the two previous cases. Here the system settles 
in a stable resonant 2:1 configuration with plan- 
ets slowly migrating inward, and eccentricities 
reaching moderate values (e ~ 0.2-0.3). 

We remark that three of the four planets 
found around the red dwarf Gliese 876 are on 
the triple 2:1 resonant configuration (Riv era et| 



al. (2010)). The masses of the planets as well as 



their distance from the star are different from 
our simulation so that our comparison is only 
qualitative but nevertheless interesting. 

To avoid spurious boundary effects, the disk 
was extended to r min = 0.5 at t — 20000, when 
the pericenter of planet #3 was close to 1. The 
radial density profile has been extrapolated from 
the inner disk edge, and the simulation restarted 
with the extended disk. 



3.4 Case Bi, 2 , 3 

In this case, cores #1 and #2 become coor- 
bital at the trap location. It is interesting that 
they remain coorbital during the growth phase 
even if they do not grow at the same time (see 
Fig. |5j. The third core is scattered out when 
planet #1 grows. It then migrates inward un- 
til it reaches the trap at a ~ 2. When its mass 
starts increasing, #3 migrates inwards, and all 
three planets reach in a compact orbital configu- 
ration. The disk has been extended to r m i„ = 0.1 
at t — 8500 using the same procedure as for case 
^3,2,1- The extended disk is followed with a time- 
step of ~ 10 -3 . To calculate 10 orbits this re- 
quires about 1 hour on 30 CPUs. The two coor- 





Q. 

< 




5000 10000 15000 20000 25000 30000 35000 40000 
1[years] 

Figure 1. Dynamical evolution of the planets in sim- 
ulation j42,i,3. Each planet is represented by differ- 
ent color: red for planet #1 (initially the innermost 
one); green for planet #2 (middle); blue for planet 
#3 (outermost). For each planet the three curves de- 
note the pericenter, semimajor axis and apocenter as 
a function of time. The masses of the cores grow to 
Jupiter mass according to Eq. |3]on At = 1000 start- 
ing t 2 (0) = 3200, ii(0) = 8000 and t 3 (0) = 12000. 
The arrows highlight t(0) for each planet. 



bital planets merge at t — 13275 and the third 
one is scattered out. The eccentricities of the two 
remaining planets grow to moderate values. The 
planets end up in the 3:1 resonance. 



3.5 Case B 3 ,i,2 

Core #3 grows first and scatters two coorbital 
cores outward. The two cores appear to be on 
the trap at very small angular separation Aa 
with: 10° < Aa < 30°. It is therefore possible 
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■8 




-8-6-4-2 8 2 4 6 8 



Figure 2. Gas surface density in the simulation 
shown in Fig. 1 at t = 15000. Three eccentric Jupiter- 
mass planets produce a complex distribution of gas. 
Color scale range was chosen such that dark blue 
corresponds to values ^ 10~ 5 and yellow to values 
> 10~ 4 (in dimensionless units). 




10000 20000 30000 40000 50000 60000 
t[years] 



Figure 3. The same as Fig.[T]for case ^3,1,2- The ar- 
rows highlight the t(0) values: t 3 (0) = 3200, ti(0) = 
8700 and t 3 (0) = 12000. 

that the scattering event affect them on the same 
way. Actually, as a result of the scattering they 
remain coorbital and only their angular separa- 
tion changes drastically. The two orbits separate 
at t = 7000 when core #1 starts growing. Once 
that happens, core #1 scatters the core #2 out- 




1 | , , , , , — 1 v - - ^ ^ 

5000 10000 15000 20000 25000 30000 35000 40000 
t[years] 



Figure 4. The same as Fig. [l] for case A3, 2,1- The 
growth of cores started at t 3 (0) = 3200, t 2 (0) = 8600 
and ii(0) = 12000, as indicated by the arrows. 



4 



B 3.5 




1 1 1 1 1 1 

5000 10000 15000 20000 

t[years] 



Figure 5. The same as Fig. Ft] for case 51,2,3- The 
growth of cores started at ii(0) = 2000, t 2 (0) = 7000 
and i 3 (0) = 9000, as indicated by the arrows. 



ward (Fig. j6j|. After the growth of core #2, a 
complex instability arises resulting in the ejec- 
tion of #3 from the system. The two remaining 
planets have moderate eccentricities and persist 
on chaotic orbits showing an outward migration 
trend till t — 33900 when the two planets merge. 
This case is comparable to case ^3,1 ,2- 

The disk has been extended to r m in — 0.1 at 
t = 8500. In order to follow the chaotic evolution 
of the two planets on reasonable CPU times, we 
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a 
< 




5000 10000 



15000 20000 
t[years] 



Figure 6. The same as Fig. [T] for case case -63,1,2- 
The growth of cores started at t3(0) = 3700, ti(0) = 
7000 and tj (0) = 9000, as indicated by the arrows. 



have then reduced the disk domain increasing 
r min to 0.3 at t = 30000. 



3.6 Summary 

From the limited number of cases that we have 
investigated so far, we can derive the following 
tentative implications for the dynamics of sys- 
tems with three giant planets and their interac- 
tion with a gas disk: 

(i) If the gas lasts long enough after the 
growth of the last giant planet, the system devel- 
ops an instability and one planet is typically lost. 
These ejections could be related to free-floating 



planets (Sumi et al. ( 2011 1) 



Indeed, only in one case out of five we obtained 
a final stable system with three giant planets. 
This is at odds with the results of Matsurnura et 
ah] ([2010 ) who found that the three-planet sys- 
tems often survive to the cud of gas disk lifetime. 
This difference may appear from the approxi- 
mate treatment of the gas disk in |Matsumura| 
et al.| (|2010|). When the gas density is strongly 



perturbed as in Fig. [2] it acts as an additional 
source of stochasticity in the planetary evolu- 
tion. Marzari et al.| ( |2010 1, who used a hydro- 
dynamical code similar to ours, also found that 
systems of fully-grown giant planets on close or- 
bits rarely survive to the end gas disk lifetime. 

(ii) The simplified two-planet systems, that 
emerge from the three-planet systems when the 
the third planet is eliminated, tend to be stabi- 
lized by their interaction with the gas disk. This 
was also pointed out in |Marzari et al.| ( |2010| ); 



Matsurnura et al. (20101 and in Moeckel, Ray 



mond, & Armitage 



(2008) 



In some cases, the 
two-planet system shows chaotic evolution till 
the system is reduced by a merging event. Future 
work on the full spatial problem will be needed to 
better explore the frequency of merging events. 

(iii) We did not find any case where the giant 
planets would end up on nearly-circular closely- 
packed orbits. This raises doubts about the ap- 
plicability of the initial conditions used in the 
models of planet scattering after the gas disk dis- 
persal (see Section 1 



(2008)) 



Chatterjee et al. (20081; 



Raymond et al. (20091; Juric and Tremaine 



25000 30000 35000 4 



GAS DISPERSAL AND GAS-FREE 
DYNAMICS 



In the previous section we assumed that the gas 
disk remains present after the accretion of the gi- 
ant planets, that is until the planets reach a sta- 
ble dynamical configuration. It is possible, how- 
ever, that the gas dispersal occurred during the 
planetary instability or soon after it, such that 
the planetary system did not have enough time 
to fully stabilize. Here we investigate these cases. 

The gas density was reduced at each time 
step dt as: 



P 



dt 



(4) 



where the coefficient r ~ 2000 years. This dis- 
sipation timescale is very short when compared 
to the 10 5 years timescale considered for photo- 



evaporation in Matsurnura et al. (20101. We find 
that, if the gas is removed too fast, planetary sys- 
tems can become unstable. In our simulations, 
we didn't observe any scattering events or merg- 
ing during the dissipation phase, so that we are 
confident that our results wouldn't change much 
using a longer dissipation timescale. 

We recall that our purpose is not to quan- 
titatively describe a specific phase of the giant- 
planet-disk interaction but to obtain a qualita- 
tive description of the whole phenomenon; this 
justify also the use of a dissipation function Q 
which is simple with respect to the description 
in Moeckel and Armitage (20121. 

When the disk gas density drops to val- 
ues below ~ 10 -10 (in dimensionless units), the 
effect of gas becomes negligible and we con- 
tinue the integration with SyMBA ( jDuncan et| 
al. (1998)). The planetary systems are evolved 
for up to 10 s years. 

We first tried a case where gas was removed 
after the planetary system has reached its final 
configuration. Fig. [7] shows the orbits for case 
^2,1,3, where the gas disk was removed in the 
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time interval [35000, 37000]. During the removal, 
planetary migration slows down and the two re- 
maining planets stay in the 2:1 resonance for 
the whole iV-body integration. The same ap- 
plies to £1,2,3, where no scattering event was 
found after the gas dispersal (gas was removed 
at [19000,20000] in this case). 

In ^3,2,1, where three planets initially sur- 
vived, they also survived for the whole length 
of the iV-body integration. The three planets re- 
mained in the 2:1 resonant chain and no scat- 
tering among them occurred (Fig. Isj> . This hap- 
pened independently of the removal time (if cho- 
sen after t = 20000 years in Fig. |4| and indepen- 
dently of the removal timescale r. 

( 2010[ ), in agreement 



In Matsumura et al. 



with our results, the two-planet systems also 
remained stable. For the three-planet systems, 

sec 



Matsumura et al. ( 2010 1 found stability (e.g. 
their Fig. 4) only in some cases. Unfortunately, 
having only one three-planet system we cannot 
test the statistical significance of our result. Note 
also that the previous versions of the SyMBA 



code used in Matsumura et al. (2010) had a later- 



identified problem when tracking closely-packed 
planetary systems (Levison et al. 2011). It has to 
be verified that this problem did not cause arti- 
ficial instabilities in some of their integrations. 

We now turn our attention to the possibility 
that the gas disk dispersed during the planetary- 
scattering phase. In Fig. [9] we removed gas in 
the interval [15000, 17000] during the evolution 
of system -4,2,1,3 (Fig. [T|. In this case, the three 
giant planets undergo a gravitational scatter- 
ing event in which one planet is ejected and 
the remaining two are sent onto highly-eccentric 
mutually-decoupled orbits (Fig. Rjb. The same 
kind of evolution happened in all cases investi- 
gated here, provided that the gas disk was re- 
moved during the scattering phase. 



5 CONCLUSIONS 

We used a hydro-dynamical code to follow the 
orbital evolution of systems of three planets as 
they grew in sequence to Jupiter mass. We found 
that the planet system changes drastically after 
the growth of each core. The orbital evolution 
of planets can be very complex. More often than 
not the orbits become unstable leading to a phase 
of planetary scattering. Planets can be ejected or 



merge ( |Marzari et al.| ( |2010[ )). 

Once the system is reduced to two planets 
the dissipative effects of gas decrease orbital ec- 
centricities of the remaining planets, and migrate 
planets into a new, stable resonant configuration. 
In only one case out of five, there was no in- 
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Figure 7. Evolution of planetary orbits in case 
Aa,i,3. The gas disk is removed at [35000, 37000], and 
an iV-body integrator is used to follow the gravita- 
tional interactions among planets for t > 37000. The 
arrow indicates the beginning of the gas-free phase. 
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Figure 8. Evolution of planetary orbits in case 
A 3 , 2 ,i. The gas disk is removed at [39500,41000]. The 
arrow indicates the beginning of the gas-free phase. 



stability happening after the growth of all three 
planets. The three-planet system remained in a 
resonant stable configuration in this case. 

If the gas disk is removed after the new sta- 
ble configuration is achieved, the orbital eccen- 
tricities remain low and the system is stable. 
This is at odds with the assumption typically 



made in the planet-scattering models (Chatter- 



jee et al.| (|2008[); |Raymond et al.| (|2009j); |Juric 

and Trcmainc (2008J)), where gas is ignored and 
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of the planetesimal disks, etc.; e.g. Tsiganis et al. 



1e+06 
1[years] 

Figure 9. Evolution of planetary orbits in case 
A 2 ,i,3. The gas disk is removed at [15000, 17000]. The 
arrow indicates the beginning of the gas-free phase. 



2005, |Levison et aL] ( |2011[ )). 

In conclusion, the results presented here 
show that the problem of understanding the dy- 
namical paths leading to the surprisingly large 
eccentricities of extrasolar planets is not fully 
resolved. Future work should improve upon our 
efforts by using more realistic prescriptions for 
the planet growth and gas dispersal, extend the 
simulation to longer timescales, and perform a 
larger number of simulations so that the statisti- 
cal significance of individual outcomes and their 
dependence on disk parameters is better under- 
stood. 
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the planets are initially placed on closely-packed 
unstable orbits. Here we show that these initial 
conditions may not naturally arise from a previ- 
ous stage, in which the planets interacted with 
their natal protoplanetary disk. 

If the gas disk is removed during the plan- 
etary instability, planetary scattering continues 
after the gas removal and the surviving planets 
can reach very eccentric orbits. However, given 
the short duration of the planetary instability 
phase, the removal of gas during this phase would 
require special timing. For example, it may be 
possible that the giant planets generally form to- 
ward the end of the disk lifetime. Or, as long 
as there is gas in the system, the existing giant 
planets keep growing and new giant planets keep 
forming. This would lead to a richer sequence of 
planetary instabilities than the one investigated 
here. We will investigate this possibility in the 
future work. 

Another possibility is that the number of gi- 
ant planets that form in a typical disk is large 
(> 3). The TV-body simulations have already 
shown that the eccentricity distribution of ex- 
oplanets implies that at least three giant planet 
existed in a typical system after the gas disk dis- 
persal. Our results seem to suggest that, for this 
condition to be fulfilled, more than three planets 
have to form originally. 

Alternatively, the giant-planet systems that 
emerge from gas disks are stable in isolation, 
as suggested by in the simulations performed in 
this work, but become unstable due to external 
causes (interactions with smaller planets, effects 
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